# Markdown
        # library(bookdown)
    # General data analysis and transformation
        library(readr)
        library(readxl)
        library(dplyr)
        library(janitor)
        # library(gtools)
        library(stringr)
        library(units)
        library(tidyr)
        library(rlang)
    # API-related
        # library(jsonlite)
        # library(urltools)
    # Mapping and GIS operations
        library(sf)
        # library(leaflet)
        # library(htmlwidgets)
        # library(geojsonsf)
        # library(rmapshaper)
        library(tmap)
        library(tmaptools)
        library(ceramic)
    # workflow
        library(here)
    # plotting
        library(ggplot2)
        library(forcats)
# Get list of CES 3 variable names and plain text names
    ces_variables <- read_csv(here('data_processed', 
                                    'ces_names.csv'))
    ces_variables <- ces_variables %>% 
        mutate(name_revised = make_clean_names(name))
    
ces_measures_overall <- ces_variables %>% 
    filter(group == 'overall', 
           subgroup == 'overall', 
           type == 'score')

ces_measures_effects <- ces_variables %>% 
    filter(group == 'pollution burden', 
           subgroup == 'environmental effects', 
           type == 'score')

ces_measures_exposures <- ces_variables %>% 
    filter(group == 'pollution burden', 
           subgroup == 'exposures', 
           type == 'score')

ces_measures_sensitive_pops <- ces_variables %>% 
    filter(group == 'population characteristics', 
           subgroup == 'sensitive populations', 
           type == 'score')

ces_measures_socioeconomic <- ces_variables %>% 
    filter(group == 'population characteristics', 
           subgroup == 'socioeconomic factors', 
           type == 'score')

Introduction

This document attempts to do some quantitative geospatial analysis to look at correlations between redline mapping and indicators of water quality, public health, and facility information. It computes the area weighed average percentile for different CalEnviroScreen 3.0 (CES) indicators by HOLC rating for a given city and statewide. CES considers two broad categories of indicators, (1) Pollution Burden and (2) Population Characteristics, which are both further broken down into two sub-categories each (factsheet, indicators page).

The Pollution Burden category includes the following indicators:

  • Environmental Effects
    • Cleanups
    • Groundwater Threats
    • Haz Waste
    • Imp Water Bodies
    • Solid Waste
  • Exposures
    • Ozone
    • PM2.5
    • Diesel PM
    • Drinking Water
    • Pesticides
    • Tox Releases
    • Traffic

The Population Characteristics category includes the following indicators:

  • Sensitive Populations
    • Asthma
    • Low Birth Weight
    • Cardiovascular Disease
  • Socioeconomic Factors
    • Education
    • Linguistic Isolation
    • Poverty
    • Unemployment
    • Housing Burden
    # Define coordinate systems to use for transformations
        projected_crs <- 3310 # see: https://epsg.io/3310 

    # load data
        # Redline Polygons
            redline_polygons <- st_read(here('data_processed', 'redline_polygons.gpkg'), quiet = TRUE)
                # st_crs(redline_polygons) # to check the reference system
            # transform to projected CRS
                redline_polygons <- redline_polygons %>% st_transform(projected_crs)
            # add the area of each redline polygon
                redline_polygons <- redline_polygons %>% mutate(holc_poly_area = st_area(.))
                
                redline_polygons <- redline_polygons %>% 
                    mutate(holc_city = factor(holc_city),
                           holc_grade = factor(holc_grade))
        # Redline cities list
            holc_cities_list <- c('Fresno', 'Los Angeles', 'Oakland',
                                  'Sacramento', 'San Diego', 'San Francisco',
                                  'San Jose', 'Stockton')
                
    # NOTE: not all redline polygons have an HOLC ID listed - create a unique ID for each one with a missing ID (so that each unique HOLC polygon can be associated with a weighted average CES score)
    # create new column for the new holc id to be assigned
        redline_polygons <- redline_polygons %>% 
            mutate('holc_id_2' = holc_id) %>% 
            arrange(holc_city, holc_grade, holc_id)
    # initialiaze some variables for the loop below
        grade_city_vec <- c()
        j <- 1
    # loop through all HOLC polygons ...
        for (i_row in 1:nrow(redline_polygons)) {
            cur_id <- redline_polygons$holc_id[i_row]
            grade <- redline_polygons$holc_grade[i_row]
            city <- redline_polygons$holc_city[i_row] 
            # the current combination of grade and city
                grade_city_cur <- paste0(city, '_', grade)
            # if the holc id is NA, check to see if there are already other missing holc id's for that city/grade combination - if so, increment a counter by 1, otherwise set counter (back) to 1
                if (is.na(cur_id)) {
                    if(grade_city_cur %in%  grade_city_vec) {
                        grade_city_counter <- grade_city_counter + 1
                    } else {
                        grade_city_counter <- 1
                    }
            # enter a new unique ID for the row with the missing ID
                redline_polygons$holc_id_2[i_row] <- paste0('unknown_', 
                                                            grade, '_', 
                                                            grade_city_counter)
            # add the grade/city combo to a running list to track number of missing holc id's for that city/grade combination
                grade_city_vec[j] <- grade_city_cur
                j <- j + 1
            }
    
    }
    
        # CES 3 Polygons
            ces3_poly <- st_read(here('data_processed', 'ces3_poly.gpkg'), quiet = TRUE)
                # st_crs(ces3_poly) # to check the reference system
            # modify the 'City' column to better keep track when joining with other datasets
                ces3_poly <- ces3_poly %>% rename('ces_city' = 'Nearby_City')
            # clean names
                ces3_poly <- ces3_poly %>% clean_names()
                
        # Transform both to projected crs
            # st_crs(redline_polygons) # redline_polygons <- redline_polygons %>% st_transform(crs = projected_crs)
            # st_crs(ces3_poly) # ces3_poly <- ces3_poly %>% st_transform(crs = projected_crs)
            
    # create a map to check that the data loaded correctly
        # tm_shape(ces3_poly) + tm_borders() + tm_shape(redline_polygons) + tm_borders(col = 'red')
    basemap <- function(city) {
    # get a basemap (See: https://github.com/mtennekes/tmap/issues/185)
        Sys.setenv(MAPBOX_API_KEY = Sys.getenv('mapbox_api_key'))
        map_background <- cc_location(ces3_poly %>% 
                                      filter(ces_city == city), 
                                  verbose = FALSE)
        map_background
    }
    holc_map <- function(city) {
        background <- basemap(city)
        # create map
        map_redline <- tm_shape(background) + 
            tm_rgb(alpha = 0.5) + 
            tm_shape(redline_polygons %>% 
                         filter(holc_city == city), is.master = TRUE) + 
            tm_polygons('holc_grade', 
                        palette = c('green', 'blue', 'yellow', 'red'), 
                        title = 'HOLC Grade') + 
            tm_layout(legend.bg.color = 'white', 
                      legend.bg.alpha = 0.7, 
                      legend.frame = TRUE)
        map_redline
    }
    # Fix self-intersecting polygons
        if (sum(!st_is_valid(ces3_poly)) > 0) {
            ces3_poly <- sf::st_buffer(ces3_poly, dist = 0)
        }
        if (sum(!st_is_valid(redline_polygons)) > 0) {
            redline_polygons <- sf::st_buffer(redline_polygons, dist = 0)
        }
    # get the intersection of the CES polygons and redline polygons
        ces_redline_intersection <- st_intersection(ces3_poly, redline_polygons)
    # inspect
        # glimpse(ces_redline_intersection)

# add the area of each clipped polygon to the data frame
        ces_redline_intersection <- ces_redline_intersection %>% 
          mutate(clipped_area = st_area(.))
            # glimpse(ces_redline_intersection)
# # write to output file 
    # st_write(obj = ces_redline_intersection,
    #          here('analysis-output', 'ces_redline_intersection.shp'))
    ces3_poly_map <- function(ces_id, ces_measure_name, city) {
        background <- basemap(city)
        # plot the CES polygons
        map_ces <- tm_shape(background) + 
            tm_rgb(alpha = 0.5) + 
            tm_shape(ces3_poly %>% 
                         filter(ces_city %in% (ces_redline_intersection %>% 
                                                   filter(holc_city == city) %>% 
                                                   pull(ces_city)))
            ) +
                         # filter(ces_city == 'Sacramento' | ces_city == 'West Sacramento')) + 
            tm_polygons(col = ces_id, 
                        title = ces_measure_name, 
                        border.alpha = 1) + 
            tm_shape(ces_redline_intersection %>% 
                         filter(holc_city == city), 
                     is.master = TRUE) +
            tm_borders(lwd = 0, alpha = 0) +
            # tm_shape(redline_polygons %>% filter(holc_city == 'Sacramento')) +
            # tm_borders(lwd = 2, col = 'black') +
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'D')) + 
            tm_borders(lwd = 2, col = 'red') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'C')) + 
            tm_borders(lwd = 2, col = 'gold2') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'B')) + 
            tm_borders(lwd = 2, col = 'blue') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'A')) + 
            tm_borders(lwd = 2, col = 'green') +
            # tm_text('holc_grade', size = 0.5) + 
            tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
        map_ces
        
    }
    ces_redline_overlap_map <- function(ces_id, ces_measure_name, city) {
        background <- basemap(city)
        map_overlap <- tm_shape(background) + 
            tm_rgb(alpha = 0.5) + 
            # tm_shape(ces3_poly %>% filter(ces_city == 'Sacramento')) + 
            tm_shape(ces3_poly %>% 
                         filter(ces_city %in% (ces_redline_intersection %>% 
                                                   filter(holc_city == city) %>% 
                                                   pull(ces_city)))
            ) +
            tm_borders(lwd = 1) + 
            tm_shape(ces_redline_intersection %>% filter(holc_city == city), is.master = TRUE) + 
            tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) + 
            # tm_shape(redline_polygons %>% filter(holc_city == 'Sacramento')) + 
            # tm_borders(lwd = 2, col = 'black') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'D')) + 
            tm_borders(lwd = 2, col = 'red') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'C')) + 
            tm_borders(lwd = 2, col = 'gold2') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'B')) + 
            tm_borders(lwd = 2, col = 'blue') + 
            tm_shape(redline_polygons %>% 
                       filter(holc_city == city, holc_grade == 'A')) + 
            tm_borders(lwd = 2, col = 'green') + 
            #tm_text('holc_grade', size = 0.5) + 
            tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
        map_overlap
    }
    # FACET MAP
    ces_redline_facet_map <- function(ces_id, ces_measure_name, city) {
        background <- basemap(city)
        map_facets <- tm_shape(ces3_poly %>% 
                         filter(ces_city %in% (ces_redline_intersection %>% 
                                                   filter(holc_city == city) %>% 
                                                   pull(ces_city)))
            ) +
            # tm_shape(ces3_poly %>% filter(ces_city == 'Sacramento')) +
            tm_borders(lwd = 1) +
            tm_shape(background) +
            tm_rgb(alpha = 0.5) +
            tm_shape(ces_redline_intersection %>% 
                         filter(holc_city == city), is.master = TRUE) +
            # tm_shape(ces_redline_intersection %>% filter(holc_city == 'Sacramento')) +
            tm_facets(by = 'holc_grade', free.coords = FALSE) +
            tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) +
            tm_shape(redline_polygons %>% filter(holc_city == city)) +
            tm_facets(by = 'holc_grade') +
            tm_borders(lwd = 2, 'black')  +
            # # D
            # tm_shape(redline_polygons %>%
            #            filter(holc_city == 'Sacramento', holc_grade == 'D')) +
            # tm_facets(by = 'holc_grade') +
            # tm_borders(lwd = 2, col = 'red') +
            # # C
            # tm_shape(redline_polygons %>%
            #            filter(holc_city == 'Sacramento', holc_grade == 'C')) +
            # tm_borders(lwd = 2, col = 'gold2') +
            # # tm_facets(by = 'holc_grade') +
            # # B
            # tm_shape(redline_polygons %>%
            #            filter(holc_city == 'Sacramento', holc_grade == 'B')) +
            # tm_borders(lwd = 2, col = 'blue') +
            # # tm_facets(by = 'holc_grade') +
            # # A
            # tm_shape(redline_polygons %>%
            #            filter(holc_city == 'Sacramento', holc_grade == 'A')) +
            # tm_borders(lwd = 2, col = 'green') +
            # tm_facets(by = 'holc_grade') +
            tm_legend(bg.color= 'white', frame = TRUE)

        map_facets
    }
ces_scores_plot <- function(ces_id, ces_measure_name) {
    
    # group by city and holc rating, then compute the weighted average score for each
        ces_redline_grouped <- ces_redline_intersection %>% 
            st_drop_geometry() %>%
            mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>% 
            group_by(holc_city, holc_grade) %>% 
            summarize(total_area = sum(clipped_area), 
                      area_x_score_total = sum(area_x_score)) %>% 
            mutate(weighted_score = area_x_score_total / total_area) %>% 
            drop_units()

    ces_redline_grouped_statewide <- ces_redline_grouped %>%
            group_by(holc_grade) %>%
            summarize(total_area = sum(total_area),
                      area_x_score_total = sum(area_x_score_total)) %>%
            mutate(weighted_score = area_x_score_total / total_area) %>%
            mutate(holc_city = '**Statewide**') %>%
            drop_units()

    ces_redline_grouped_city_state <- bind_rows(ces_redline_grouped, ces_redline_grouped_statewide)

    # plot the results by city and holc rating, and include statewide summary
        g_city_state <- ggplot(ces_redline_grouped_city_state) +
            aes(x = fct_relevel(fct_reorder(holc_city, weighted_score), 
                                c('**Statewide**'), 
                                after = 0), 
                y = weighted_score, 
                fill = fct_rev(holc_grade)) +
            geom_bar(stat = 'identity', position = 'dodge') +
            geom_vline(xintercept = 1.52, size = 0.5, color = 'grey50') +
            scale_fill_manual(values = rev(c('green', 'blue', 'yellow', 'red'))) +
            guides(fill = guide_legend(reverse = TRUE)) +
            labs(fill = 'HOLC Rating', 
                 title = ces_measure_name,
                 x = NULL, 
                 y = "Area Weighted Average Score By City") +
            coord_flip()

        g_city_state
}
ces_summary_plots <- function(ces_id, ces_measure_name) {
    # # TROUBLESHOOTING
    #     ces_measures_effects_scores <- ces_variables %>%
    #         filter(group == 'pollution burden',
    #                subgroup == 'environmental effects',
    #                type == 'score')
    #     i_map <- 4
    #     ces_measure_name <- ces_measures_effects_scores$name[i_map]
    #     ces_id <- ces_measures_effects_scores$name_revised[i_map]
    
    
    var_name <- ces_id # paste0('weighted_score_', ces_id)
    
    # compute weighted average score for each HOLC ID (first compute for ces_redline_intersection polygon
        holc_poly_wt_score <- ces_redline_intersection %>% 
            st_drop_geometry() %>%
            mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>% # compute weighted avg score for each ces_redline_intersection polygon
            group_by(holc_city, holc_id_2, holc_grade) %>% # compute weighted avg score for each HOLC polygon (combining individual ces_redline_intersection polygons that are within each HOLC polygon)
            summarize(total_area = sum(clipped_area), 
                      area_x_score_total = sum(area_x_score), 
                      component_poly_count = n()) %>% 
            mutate(!!var_name := area_x_score_total / total_area) %>% 
            ungroup() %>% 
            drop_units()
            
    # for each HOLC polygon, compute departure of weighted average score from the citywide average for the city it's in
    holc_poly_wt_score_departures <- holc_poly_wt_score %>% 
        group_by(holc_city) %>% 
        mutate(!!paste0(var_name, '_city_avg') := mean(!!as.name(var_name))) %>%
        # summarize(city_avg = mean(!!as.name(var_name))) %>% 
        ungroup() %>%
        mutate(!!paste0(var_name, '_departure') := !!as.name(var_name) - !!as.name(paste0(var_name, '_city_avg'))) %>% 
        select(-total_area, -area_x_score_total, -component_poly_count) %>% 
        {.}
    
    # compute z scores for each holc polygon, by city
    holc_z_scores <- holc_poly_wt_score_departures %>% 
        group_by(holc_city) %>% 
        mutate(!!paste0(var_name, '_sd_city') := sd(!!as.name(paste0(var_name, '_departure')))) %>% 
        mutate(!!paste0(var_name, '_z_score_city') := !!as.name(paste0(var_name, '_departure')) / !!as.name(paste0(var_name, '_sd_city'))) %>% 
        ungroup() %>% 
        {.}
    
    # # check - summary of average and median z scores by city
    #     holc_z_scores %>% 
    #         group_by(holc_city) %>% 
    #         summarise(med = round(median(cleanups_z_score_city), digits = 3), 
    #                   avg = round(mean(cleanups_z_score_city), digits = 3))
    
    
        # # DELETE
        # # box plot of z scores
        #     delete_boxplot <- ggplot(holc_z_scores, 
        #                                  aes(x = holc_city, 
        #                                      y = !!as.name(paste0(var_name, '_z_score_city')))) + 
        #         geom_boxplot()
        #     delete_boxplot
    
    
    # # for info purposes, compute the statewide average of each holc grade's departure
    #     statewide_avg_departures_bygrade <- holc_poly_wt_score_departures %>% 
    #         group_by(holc_grade) %>% 
    #         summarize(sw_avg_dep = mean(!!as.name(paste0(var_name, '_departure'))),
    #                   sw_median_dep = mean(!!as.name(paste0(var_name, '_departure'))))

    # box plot - raw scores
    raw_scores_boxplot <- ggplot(data = holc_poly_wt_score_departures, 
                                 aes(x = holc_grade, 
                                     y = !!as.name(paste0(var_name)))) + 
        geom_boxplot(aes(fill = holc_grade), notch = TRUE) + 
        # scale_color_manual(values = alpha(c('green', 'blue', 'yellow', 'red'), 0.4)) +
        scale_fill_manual(values = alpha(c('green', 'blue', 'yellow', 'red'), 0.6)) + 
        geom_jitter(color='black', size=0.6, alpha=0.5, width = 0.2) + 
        theme(legend.position="none") + 
        labs(x = 'HOLC Grade', y = paste0(ces_measure_name, ' Raw Score')) + 
        geom_blank()
    # raw_scores_boxplot

    # box plot - departures
        departures_boxplot <- ggplot(data = holc_poly_wt_score_departures, 
                                     aes(x = holc_grade, 
                                         y = !!as.name(paste0(var_name, '_departure')))) + 
            geom_boxplot(aes(fill = holc_grade), notch = TRUE) + 
            # scale_color_manual(values = alpha(c('green', 'blue', 'yellow', 'red'), 0.4)) +
            scale_fill_manual(values = alpha(c('green', 'blue', 'yellow', 'red'), 0.6)) + 
            geom_jitter(color='black', size=0.6, alpha=0.5, width = 0.2) + 
            theme(legend.position="none") + 
            labs(x = 'HOLC Grade', y = paste0(ces_measure_name, ' Departure')) + 
            geom_blank()
        # departures_boxplot

    # # density plot
    #     departures_density <- ggplot(holc_poly_wt_score_departures, 
    #                                  aes(x = !!as.name(paste0(var_name, '_departure')))) + 
    #         geom_density()
    #     departures_density
    
    # point plot - all cities - departures
        departures_point <- ggplot(holc_poly_wt_score_departures) + 
            geom_jitter(aes(x = !!as.name(paste0(var_name, '_departure')), 
                           y = holc_city, color = holc_grade), 
                        height = 0.25) + 
            scale_color_manual(values = alpha(c('green', 'blue', 'orange', 'red'), 0.3)) +
            labs(x = paste0(ces_measure_name, ' Departure'), y = 'City') +
            scale_y_discrete(limits = rev(levels(holc_poly_wt_score_departures$holc_city))) + 
            geom_blank()
        # departures_point
        
    # point plot - all cities - raw scores
        raw_scores_point <- ggplot(holc_poly_wt_score_departures) + 
            geom_jitter(aes(x = !!as.name(var_name), 
                           y = holc_city, color = holc_grade), 
                        height = 0.25) + 
            scale_color_manual(values = alpha(c('green', 'blue', 'orange', 'red'), 0.3)) +
            labs(x = paste0(ces_measure_name, ' Raw Score'), y = 'City') +
            scale_y_discrete(limits = rev(levels(holc_poly_wt_score_departures$holc_city))) + 
            geom_blank()
        # raw_scores_point
        
    # point plot - departures - by city (facet on cities)
        departures_point_city <- ggplot(holc_poly_wt_score_departures %>% 
                                            #filter(holc_city == 'Sacramento') %>% 
                                            {.}) + 
            geom_point(aes(x = !!as.name(paste0(var_name, '_departure')), 
                           y = holc_grade, 
                           color = holc_grade)) + 
            scale_color_manual(values = alpha(c('green', 'blue', 'orange', 'red'), 0.4)) +
            scale_y_discrete(limits = rev(levels(holc_poly_wt_score_departures$holc_grade))) +
            facet_wrap(vars(holc_city)) + 
            labs(x = paste0(ces_measure_name, ' Departure'), y = 'HOLC Grade') +
            geom_blank()
        # departures_point_city
        
    # point plot - raw scores - by city (facet on cities)
        raw_scores_point_city <- ggplot(holc_poly_wt_score_departures %>% 
                                            #filter(holc_city == 'Sacramento') %>% 
                                            {.}) + 
            geom_point(aes(x = !!as.name(paste0(var_name)), 
                           y = holc_grade, 
                           color = holc_grade)) + 
            scale_color_manual(values = alpha(c('green', 'blue', 'orange', 'red'), 0.4)) +
            scale_y_discrete(limits = rev(levels(holc_poly_wt_score_departures$holc_grade))) +
            facet_wrap(vars(holc_city)) + 
            labs(x = paste0(ces_measure_name, ' Raw Score'), y = 'HOLC Grade') +
            geom_blank()
        # raw_scores_point_city
        
# compute average departure for each HOLC grade within each city
    holc_poly_wt_score_departures_summary <- holc_poly_wt_score_departures %>% 
        group_by(holc_city, holc_grade) %>% 
        summarize(city_holc_dep_average = mean(!!as.name(paste0(var_name, '_departure'))),
                  city_holc_dep_median = median(!!as.name(paste0(var_name, '_departure'))))
    
    holc_city_means_plot <- ggplot(holc_poly_wt_score_departures_summary) + 
        geom_point(aes(x = city_holc_dep_average, 
                       y = holc_city, 
                       color = holc_grade)) + 
        # geom_jitter(aes(x = city_holc_dep_average, 
        #                 y = holc_city, 
        #                 color = holc_grade),
        #             height = 0.15) + 
        scale_color_manual(values = alpha(c('green', 'blue', 'orange', 'red'), 1.0), 
                           name = 'HOLC Grade') +
        scale_y_discrete(limits = rev(levels(holc_poly_wt_score_departures_summary$holc_city))) +
        labs(x = paste0(ces_measure_name, ' Average Departure'), y = 'City') +
        theme(legend.position = 'bottom') +
        geom_blank()
    # holc_city_means_plot
    
    holc_city_medians_plot <- ggplot(holc_poly_wt_score_departures_summary) +
        geom_point(aes(x = city_holc_dep_median, y = holc_city, color = holc_grade))
    # holc_city_medians_plot
    # return the plots from the function
        list(departures_boxplot, departures_point, raw_scores_point, 
             departures_point_city, raw_scores_point_city, holc_city_means_plot,
             raw_scores_boxplot)
}

CalEnviroScreen Overall Scores

CES 3 Score

Maps

Fresno
    i_city <- 1
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

Los Angeles
    i_city <- 2
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

Oakland
    i_city <- 3
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

Sacramento
    i_city <- 4
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

San Diego
    i_city <- 5
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

San Francisco
    i_city <- 6
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

San Jose
    i_city <- 7
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

Stockton
    i_city <- 8
    ces_measure_name <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name) %>%
        {.}
    ces_id <- ces_variables %>%
        filter(id == 'CIscore') %>%
        pull(name_revised) %>%
        {.}

    map_redline <- holc_map(city = holc_cities_list[i_city])
    map_ces <- ces3_poly_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name, city = holc_cities_list[i_city])

Plots

    ces_measure_name <- ces_variables %>% 
        filter(id == 'CIscore') %>% 
        pull(name) %>% 
        {.}
    ces_id <- ces_variables %>% 
        filter(id == 'CIscore') %>% 
        pull(name_revised) %>% 
        {.}

    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    
    # box plots
    summary_plots[[7]]

    summary_plots[[1]]

    # point plots by city
    summary_plots[[3]]

    summary_plots[[2]]

    # facet plots
    summary_plots[[5]]

    summary_plots[[4]]

    # average scores by grade
    summary_plots[[6]]

Pollution Burden - Environmental Effects

Cleanups

    i_map <- 1
    ces_measure_name <- ces_measures_effects$name[i_map]
    ces_id <- ces_measures_effects$name_revised[i_map]

    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Groundwater Threats

    i_map <- 2
    ces_measure_name <- ces_measures_effects$name[i_map]
    ces_id <- ces_measures_effects$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)

    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Haz Waste

    i_map <- 3
    ces_measure_name <- ces_measures_effects$name[i_map]
    ces_id <- ces_measures_effects$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Imp Water Bodies

    i_map <- 4
    ces_measure_name <- ces_measures_effects$name[i_map]
    ces_id <- ces_measures_effects$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Solid Waste

    i_map <- 5
    ces_measure_name <- ces_measures_effects$name[i_map]
    ces_id <- ces_measures_effects$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)

    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Pollution Burden - Exposures

Ozone

    i_map <- 1
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

PM2.5

    i_map <- 2
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Diesel PM

    i_map <- 3
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Drinking Water

    i_map <- 4
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Pesticides

    i_map <- 5
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Tox Releases

    i_map <- 6
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Traffic

    i_map <- 7
    ces_measure_name <- ces_measures_exposures$name[i_map]
    ces_id <- ces_measures_exposures$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Population Characteristics - Sensitive Populations

Asthma

    i_map <- 1
    ces_measure_name <- ces_measures_sensitive_pops$name[i_map]
    ces_id <- ces_measures_sensitive_pops$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Low Birth Weight

    i_map <- 2
    ces_measure_name <- ces_measures_sensitive_pops$name[i_map]
    ces_id <- ces_measures_sensitive_pops$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Cardiovascular Disease

    i_map <- 3
    ces_measure_name <- ces_measures_sensitive_pops$name[i_map]
    ces_id <- ces_measures_sensitive_pops$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Population Characteristics - Socioeconomic Factors

Education

    i_map <- 1
    ces_measure_name <- ces_measures_socioeconomic$name[i_map]
    ces_id <- ces_measures_socioeconomic$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Linguistic Isolation

    i_map <- 2
    ces_measure_name <- ces_measures_socioeconomic$name[i_map]
    ces_id <- ces_measures_socioeconomic$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Poverty

    i_map <- 3
    ces_measure_name <- ces_measures_socioeconomic$name[i_map]
    ces_id <- ces_measures_socioeconomic$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Unemployment

    i_map <- 4
    ces_measure_name <- ces_measures_socioeconomic$name[i_map]
    ces_id <- ces_measures_socioeconomic$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

Housing Burden

    i_map <- 5
    ces_measure_name <- ces_measures_socioeconomic$name[i_map]
    ces_id <- ces_measures_socioeconomic$name_revised[i_map]
    
    # map_ces <- ces3_poly_map(ces_id, ces_measure_name)
    # map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    # 
    # tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
    # 
    # ces_redline_facet_map(ces_id, ces_measure_name)
    # 
    # ces_scores_plot(ces_id, ces_measure_name)
    
    summary_plots <- ces_summary_plots(ces_id, ces_measure_name)
    summary_plots[[1]]

    summary_plots[[4]]

# ces_measures_combined <- bind_rows(ces_measures_overall,
#                                    ces_measures_effects, 
#                                    ces_measures_exposures, 
#                                    ces_measures_sensitive_pops, 
#                                    ces_measures_socioeconomic)
# 
# for (i_output in seq(nrow(ces_measures_combined))) {
#     ces_measure_name <- ces_measures_combined$name[i_output]
#     ces_id <- ces_measures_combined$name_revised[i_output]
#     var_name <- ces_id # paste0('weighted_score_', ces_id)
#     
#     # group by city and holc ID, then compute the weighted average score for each holc ID
#         ces_redline_grouped_holc_id <- ces_redline_intersection %>% 
#             mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>% 
#             group_by(holc_city, holc_id_2, holc_grade) %>% 
#             summarize(total_area = sum(clipped_area), 
#                       area_x_score_total = sum(area_x_score), 
#                       component_poly_count = n()) %>% 
#             mutate(!!var_name := area_x_score_total / total_area) %>% 
#             drop_units()
#         
#         ces_redline_grouped_holc_id <- ces_redline_grouped_holc_id %>% 
#             select(holc_city, holc_id_2, holc_grade, !!var_name) %>% 
#             rename(city = holc_city)
#     
#     # add to output df
#         if (i_output == 1) {
#             df_output <- ces_redline_grouped_holc_id
#         } else {
#             df_output <- left_join(x = df_output, 
#                                    y = ces_redline_grouped_holc_id %>% 
#                                        st_drop_geometry(), 
#                                    by = c('city', 'holc_id_2', 'holc_grade'))
#         }
# } 
# 
# df_output <- df_output %>% arrange(city, holc_grade, holc_id_2)
# 
# 
# # # write the output shapefile
# #     st_write(obj = df_output, 
# #              here('analysis-output', 
# #                   'redline_ces_scores',
# #                   paste0('redline_ces_scores', 
# #                                             '_', 
# #                                             Sys.Date(),
# #                                             '.shp'))
# #     )
#     
# # write the output to a geopackage file
#     st_write(obj = df_output, 
#              here('analysis-output', 
#                   'redline_ces_scores',
#                   paste0('redline_ces_scores', 
#                                             '_', 
#                                             Sys.Date(),
#                                             '.gpkg'))
#     )